Непараметрические модели прогнозирования

Экспоненциальное сглаживание

Пусть \(y_t\) одномерный временной ряд. Экспоненциальным сглаживанием на период \(t+h\), основываясь на знании значений процесса до момента времени \(t\),называется построение прогноза по переменной \(\widehat a_t\). \[\widehat{y}_{t+h}=\widehat a_t\],
Переменная \(\widehat a_t\) рекурсивно пересчитывается. \[\widehat a_t= \alpha y_t+(1-\alpha)\widehat y_t=\alpha y_t+(1-\alpha)\widehat a_{t-1}\] где \(0 < \alpha <1\) называется параметром сглаживания. Начальное значение \(\widehat a_1= y_1\). Эквивалентная запись модели прогнозирования называется error correction form: \[\widehat a_t= \widehat a_{t-1}+\alpha (y_t-\widehat a_{t-1})=\widehat a_{t-1}+ \alpha (y_t- \widehat y_t) = \widehat a_{t-1} +\alpha \epsilon_t\] Экспоненциальное сглаживание это частный случай метода Хольта-Уинтерса (Holt-Winters) для несезонных рядов без ярко выраженного тренда. Посмотрим пример прогнозирование ряда.

promData <- read.csv("c:/ShurD/aLection/PromProizv.csv")
promData <- ts(promData,start= c(1995,1),frequency = 12)
plot.ts(promData,col = "blue",lwd = 2,type = "l",main = "Prom.Proizvod.")

past <- window(promData,end = c(2005,1))
future <- window(promData,start = c(2005,2))
alpha <-0.2
model <- HoltWinters(past,alpha = alpha,beta = FALSE, gamma= FALSE)
filter(alpha*past,filter = 1-alpha,method = "recursive",init = past[1])
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 1995  99.50000 100.74000  99.97200  99.53760  99.65008  99.48006 100.02405
## 1996  99.95515 100.86412 100.77130  98.91704  97.29363  96.35490  95.46392
## 1997  96.79765  97.93812  98.21049  96.90840  95.66672  95.65337  95.96270
## 1998  99.67110 100.75688 100.66550  98.35240  96.28192  94.42554  92.72043
## 1999  93.33148  94.94518  96.25615  98.22492 100.37993 102.86395 105.49116
## 2000 110.80084 110.56067 109.54854 109.75883 109.76706 109.51365 109.65092
## 2001 105.96000 105.48800 105.43040 105.74432 105.33546 105.16836 105.15469
## 2002 103.58418 103.60734 103.74588 103.55670 103.72536 104.54029 104.31223
## 2003 104.38588 104.84871 105.29896 105.93917 106.15134 106.34107 106.17286
## 2004 107.41053 107.24842 107.13874 106.81099 107.28879 106.71103 106.72883
## 2005 104.86482                                                            
##            Aug       Sep       Oct       Nov       Dec
## 1995 100.27924 101.20339 101.18271 101.98617 100.66894
## 1996  95.27114  96.51691  96.27353  97.65882  96.84706
## 1997  96.43016  98.28413  99.18730 101.12984 100.16387
## 1998  91.27634  90.80108  90.82086  91.33669  92.58935
## 1999 108.43293 108.80634 109.62507 109.92006 110.07605
## 2000 109.16074 109.40859 109.04687 107.73750 107.25000
## 2001 104.88375 104.92700 104.88160 104.42528 103.98023
## 2002 104.54978 104.41983 103.69586 103.59669 103.85735
## 2003 106.53828 106.67063 106.75650 106.98520 107.08816
## 2004 106.08306 105.56645 105.65316 105.48253 104.80602
## 2005
pred <- predict(model,n.ahead =100)
plot(model,predicted.values = pred,lwd=2)
lines(future,col = "blue")

print(model)
## Holt-Winters exponential smoothing without trend and without seasonal component.
## 
## Call:
## HoltWinters(x = past, alpha = alpha, beta = FALSE, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 0.2
##  beta : FALSE
##  gamma: FALSE
## 
## Coefficients:
##       [,1]
## a 104.8648

Мы построили прогноз на 100 месяцев вперед. Понятно, что прогноз будет константой, так как модель не подразумевает тренда.

Модель Хольта-Уинтерса с линейным трендом

Очевидным обобщением метода является добавление тренда, начнем с линейного. Модель для прогнозирования будет теперь такой \[\widehat{y}_{t+h}=\widehat a_t+h\widehat b_t\], где пересчет параметров \(\widehat a_t,\widehat b_t\) идет уже по следующим формулам \[\widehat a_t = \alpha y_t+(1-\alpha)(\widehat a_{t-1}+\widehat b_{t-1})\] \[\widehat b_t = \beta (\widehat a_t- \widehat a_{t-1})+(1-\beta)\widehat b_{t-1}\]

past <- window(promData,end = c(2011,3))
future <- window(promData,start = c(2011,4))
model <- HoltWinters(past, gamma= FALSE)
pred <- predict(model,n.ahead =40)
plot(model,predicted.values = pred,lwd=2)
lines(future,col = "blue",lwd = 3)

print(model)
## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = past, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 0.7866934
##  beta : 0.1602257
##  gamma: FALSE
## 
## Coefficients:
##          [,1]
## a 104.6967136
## b  -0.1611357

Коэффициенты \(\alpha ,\beta\) мы не задавали. Они найдены были оптимизацией в смысле минимума квадрата ошибки прогноза на один шаг вперед по первой части ряда.

Cезонный метод Хольта-Уинтерса

Дальнейшее расширение метода это добавление аддитивной или мультипликативной сезонной компоненты с периодом сезонности \(p\).

Аддитивная сезонная компонента

Модель будет теперь следующего вида. \[\widehat{y}_{t+h}=\widehat a_t+h\widehat b_t+ \widehat s_t\], где пересчет параметров \(\widehat a_t,\widehat b_t,\widehat s_t\) идет уже по следующим формулам \[\widehat a_t = \alpha (y_t-\widehat s_{t-p})+(1-\alpha)(\widehat a_{t-1}+\widehat b_{t-1})\] \[\widehat b_t = \beta (\widehat a_t- \widehat a_{t-1})+(1-\beta)\widehat b_{t-1}\] \[\widehat s_t = \gamma (y_t-\widehat a_t)+(1-\gamma)\widehat s_{t-p}\]

Мультипликативная сезонная компонента

Модель будет теперь следующего вида. \[\widehat{y}_{t+h}=(\widehat a_t+h\widehat b_t)\widehat s_t\], где пересчет параметров \(\widehat a_t,\widehat b_t,\widehat s_t\) идет уже по следующим формулам \[\widehat a_t = \alpha (y_t-\widehat s_{t-p})+(1-\alpha)(\widehat a_{t-1}+\widehat b_{t-1})\] \[\widehat b_t = \beta (\widehat a_t- \widehat a_{t-1})+(1-\beta)\widehat b_{t-1}\] \[\widehat s_t = \gamma (y_t/\widehat a_t)+(1-\gamma)\widehat s_{t-p}\]

Числовой пример рассмотрим на примере ряда импорта России из прошлой лекции

import <- read.csv("c:/ShurD/Alection/import.csv", header = TRUE)
importFact <- ts(import$Fact ,start= c(1997,1),frequency =12)
plot(importFact ,type = "l", col = "blue",lwd = 2,main = "Import (mlrd dollar).Russia")

Начнем прогнозирование с января 2011 года на 5 лет вперед или на 12*5 = 60 месяцев вперед. Модель будем брать мультипликативную.

past <- window(importFact,end = c(2011,1))
future <- window(importFact,start = c(2011,2))
model <- HoltWinters(past,seasonal = "mult")
pred <- predict(model,n.ahead =60)
plot(model,predicted.values = pred,lwd=2)
lines(future,col = "blue",lwd = 3)

print(model)
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = past, seasonal = "mult")
## 
## Smoothing parameters:
##  alpha: 0.6389659
##  beta : 0.01602396
##  gamma: 0.8568556
## 
## Coefficients:
##           [,1]
## a   24.4933554
## b    0.1761019
## s1   0.8849390
## s2   1.0279983
## s3   1.0611034
## s4   1.0875641
## s5   1.1375211
## s6   1.1930511
## s7   1.2343608
## s8   1.1516480
## s9   1.1371452
## s10  1.0515990
## s11  1.1456385
## s12  0.6530003

2012 год спрогнозировали неплохо, а потом запахло новым кризисом:нефть стала падать,в 2015 году Украина, “Крым наш”, санкции ввели, нефть совсем упала и т.д.. Вывод. Прогнозировать на большой период имеет смысл только при неизменной динамике ряда.

Помочь увидеть недостатки и понять область применимости метода Хольта-Уинтерса поможет следующий фрейм

Параметрические модели прогнозирования

Чаще всего построение параметрической модели временных рядов вызвано желанием успешного прогнозирования будущих значений.В этом разделе будем предполагать, что модель уже известна.

Прогноз с минимальной среднеквадратичной ошибкой

Через \(\hat{y}_t(l)\) будем обозначать прогноз,сделанный в момент времени \(t\) на \(l\) интервалов времени вперед, а \(y_{t+l}\) истинное значение временного ряда, которое он примет в будущем через \(l\) интервалов времени после \(t\).

\(\hat y_{t}\) можно выразить как линейную функцию от предшествующих значений белого шума \(\epsilon_t,\epsilon_{t-1},....\)

Пусть наилучшим является прогноз вида \(\hat y_{t}(l)=\Psi_{l}^{∗}\epsilon_{t}+\Psi_{l+1}^{∗}\epsilon_{t-1}+\Psi_{l+2}^{∗}\epsilon_{t-2}+...=\sum_{j=0}^{∞}\Psi_{l+j}^{∗}\epsilon_{t-j}\)

где веса \(\Psi_{l}^{∗},\Psi_{l+1}^{∗},\Psi_{l+2}^{∗},..\) пока неизвестны их необходимо определить. Тогда так как

\(y_{t+l}=\sum_{j=-∞}^{t+l}\Psi_{t+l-j}\epsilon_{j}=\sum_{j=0}^{∞}\Psi_{j}\epsilon_{t+l-j}\)

Вычислим математическое ожидание квадрата ошибки прогноза

\(E[y_{t+l}-\hat y_{t}(l)]^2=\sigma_{\epsilon}^2(1+\Psi_1^2+...+\Psi_{l-1}^2)+\sigma_{\epsilon}^2\sum_{j=0}^{∞}(\Psi_{l+j}-\Psi_{l+j}^{∗})^2\)

Оба слагаемых в правой части последнего выражения положительны. Второе слагаемое можно сделать нулем если положить \(\Psi_{l+j}=\Psi_{l+j}^{∗}\) Таким образом при \(\Psi_{l+j}=\Psi_{l+j}^{∗}\) математиматическое ожидание квадрата ошибки минимально. Так мы доказали частный случай важного свойства условного математического ожидания

Условное математическое ожидание \(\mathbb{E}[X \mid \mathcal{G}]\) — это наилучшее средне-квадратичное приближение \(X \mathcal{G}\)-измеримыми случайными величинами:

\(\|X - \mathbb{E}[X\mid \mathcal{G}]\| = \inf\limits_{Z \in L^2_{\mathcal{G}}} \|X - Z\|\)

См. свойства условного математического ожидания например здесь https://ru.wikipedia.org/wiki/%D0%A3%D1%81%D0%BB%D0%BE%D0%B2%D0%BD%D0%BE%D0%B5_%D0%BC%D0%B0%D1%82%D0%B5%D0%BC%D0%B0%D1%82%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%BE%D0%B5_%D0%BE%D0%B6%D0%B8%D0%B4%D0%B0%D0%BD%D0%B8%D0%B5 Итак

\(\hat y_t(l)=E[y_{t+l}|y_1,..y_t]\) - наилучший в смысле среднего квадрата ошибки прогноз. Вычислением и изучением свойст прогноза мы и будем заниматься в этой лекции для различных моделей.

Детерминированный тренд

Пусть модель временного ряда

\(y_t=\mu_t+x_t\)

где случайная компонента \(x_t\) имеет среднее 0, а \(\mu_t\)- детерминированный тренд

\(\hat y_t(l)=E[\mu_{t+l}+X_{t+l}|y_1,..y_t]=E[\mu_{t+l}|y_1,..y_t]+E[x_{t+l}|y_1,..y_t]=\mu_{t+l}+E[x_{t+l}]=\mu_{t+l}\)

так как при \(l>1\) \(x_{t+l}\) не зависит от \(y_1,..y_t\) и имеет среднее 0

В частности для линейного тренда \(\mu_t=\beta_0+\beta_1t\) имеем прогноз

\(\hat y_t(l)= \beta_0+\beta_1(t+l)\)

Для сезонной модели \(\mu_t=\mu_{t+12}\)

\(\hat y_t(l)= \mu_{t+12+l}\)

Ошибка прогноза \(e_t(l)= y_{t+l}-\hat y_t(l)= \mu_{t+l}+x_{t+l}-\mu_{t+l}= x_{t+l}\)

Отсюда \(E[e_t(l)]=E[x_{t+l}]= 0\) и \(D[e_t(l)]= c_0 = D[x_t]\)

ARIMA прогнозирование

AR(1)

На примере AR(1) модели с константой

\(y_t-\mu=\phi(y_{t-1}-\mu)+\epsilon_t\)

Заменим \(t\) на \(t+1\)

\(y_{t+1}-\mu=\phi(y_{t}-\mu)+\epsilon_{t+1}\)

возьмем условное матеиатическое ожидание от обеих частей

\(\hat y_t(1)-\mu = \phi[E(y_t|y_1,...,y_t)-\mu]+[E(\epsilon_{t+1}|y_1,...,y_t)\)

Отсюда

\(\hat y_t(1)=\mu+\phi(y_t-\mu)\)

Для произвольного \(l\)

\(\hat y_t(l)=\mu+\phi(\hat y_t(l-1)-\mu)\) для \(l >=1\)

Последнее можно переписать так

\(\hat y_t(l)=\mu+\phi^{l-1}(\hat y_t(1)-\mu)\) для \(l >=1\)

или

\(\hat y_t(l)=\mu+\phi^{l}(y_t-\mu)\)

Посмотрим ошибку при прогнозировании

\(e_t(1)=y_{t+1}-\hat{y}_t(1)= \phi(y_t-\mu)+\mu+\epsilon_{t+1}-\phi(y_t-\mu)+\mu= \epsilon_{t+1}\)

Отсюда дисперсия ошибки прогноза

\(D[e_t(1)]=\sigma_{\epsilon}^2\)

Ошибка прогноза на \(l\) интервалов времени

\(e_t(l)=y_{t+l}-\hat{y}_t(l1)=\epsilon_{t+l}+\phi\epsilon_{t+l-1}+...+\phi^{l-1}\epsilon_{t+1}\)

Но для модели AR(1) \(\Psi_k= \phi^k\) поэтому

\(e_t(l)=\epsilon_{t+l}+\Psi_1\epsilon_{t+l-1}+...+\Psi_{l-1}\epsilon_{t+1}\)

следовательно дисперсия ошибки прогноза на \(l\) интервалов

\(D[e_t(l)]= \sigma_{\epsilon}^2(1+\Psi_1^2+...+\Psi_{l-1}^2)=\sigma_{\epsilon}^2\frac{1-\phi^{2l}}{1-\phi^2}\)

MA(1)

Модель

\(y_t=\mu+\epsilon_t-\theta\epsilon_{t-1}\)

Меняем \(t\) на \(t+1\) и ,беря условное математическое ожидание от обеих часте, получим

\(\hat{y}_t(1)= \mu-\theta E[\epsilon_t|y_1,...,y_t]=\mu-\theta\epsilon_t\)

Откуда ошибка прогноза

\(e_t(1)=y_{t+1}-\hat{y}_t(1)= (\mu+\epsilon_{t+1}-\theta\epsilon_{t})- (\mu-\theta\epsilon_t)= \epsilon_{t+1}\)

Для произвольного числа прогнозов

\(\hat{y}_t(l)= \mu+ E[\epsilon_{t+l}|y_1,...,y_t]- \theta E(\epsilon_{t+l-1}|y_1,...,y_t)=\mu\) для всех\(l>1\)

Для MA(1) модели \(\Psi_1=-\theta\) и \(\Psi_j=0\) при \(j>1\) и формула

\(D[e_t(l)]= \sigma_{\epsilon}^2(1+\Psi_1^2+...+\Psi_{l-1}^2)\)

также справедлива.

ARMA(p,q)

Для модели произвольного порядка уравнение для прогнозирования

\(\hat{y}_t(l)= \phi_1\hat{y}_{t-1}+\phi_2\hat{y}_{t-2}+...+\phi_p\hat{y}_{t-p}+\theta_0-\theta_1E(\epsilon_{t+l-1}|y_1,...,y_t)-\theta_2E(\epsilon_{t+l-2}|y_1,...,y_t)-...--\theta_qE(\epsilon_{t+l-q}|y_1,...,y_t)\)

где

\(E(\epsilon_{t+j}|y_1,...,y_t)= 0\) при \(j>0\)

\(E(\epsilon_{t+j}|y_1,...,y_t)= \epsilon_{t+j}\) при \(j\le 0\)

Заметим, что для \(j>0\) это будет прогноз \(\y_{t}(j)\), а для \(j\le 0\) это известное уже значение \(\hat{y}_(j)= y_{t+j}\)

Немного сложнее, но также можно показать, что дисперсии ошибок справедливо соотношение

\(D[e_t(l)]= \sigma_{\epsilon}^2(1+\Psi_1^2+...+\Psi_{l-1}^2)\)

Дверительный интервал прогноза

Здесь и понадобится выражение для дисперсии ошибок, полученное в предыдущих разделах

\(D[e_t(l)]= \sigma_{\epsilon}^2(1+\Psi_1^2+...+\Psi_{l-1}^2)\) (*)

Типичным предположением относительно распределения белого шума это нормальность. На практике \(\sigma_\epsilon^2\) неизвестна и должна быть оцененена из рассматриваемого временного ряда. Пусть порядки ARMA модели \(p\) и \(q\) определены правильно.Предположим, что ряд стационарен. Одним из рассмотренных ранее методом произведем оценку параметров модели \[\phi_1,...\phi_p,\theta_1,...,\theta_q\] и \(\sigma_\epsilon^2\). Для нахождения доверительного интервала с использованием формулы для дисперсии необходимо определить веса \(\Psi_k, k=1,...l-1\). Это мы ранее сделали, для моделей AR(1) и МА(1) в явном виде. Для произвольной модели ARMA(p,q) это можно сделать численно. Имеем модель \[\phi_p(B)y_t = \theta_q(B)\epsilon_t\]. Где \[\phi_p(B)= 1- \phi_1 B- ...- \phi_p B^p\] , а \[\theta_q(B)=1-\theta_1 B-...-\theta_qB^q\] Реализовав численный алгоритм деления полиномов \[\frac{\phi_p(B)}{\theta_q(B)}\], найдем необходимые веса \(\Psi_k, k=1,...l-1\). Таким образом дисперсия прогноза на \(l\) интервалов нам известна это формула (*). Для заданного уровня доверия \((1-\alpha)100\)% определим квантиль \(z_{1-\alpha/2}\) стандартного нормально распределения с нулевым математическим ожиданием и дисперсией 1, тогда доверительный интервал для прогноза \(\hat y_t(l)\) будет

\(\hat y_t(l)\pm z_{1-\alpha/2}\sqrt{D[e_t(l)]}\)

Пример cмоделированного ряда прогнозирования

Смоделируем ARMA(2,1) модель

\((1-\phi_1B-\Phi_2B^2)y_t= (1-\theta_1B)\epsilon_t\)

phi1 <- 0.4
phi2 <- -0.5
theta1 <- 0.4
sigma <- 4
model <- list(ar = c(phi1,phi2) ,ma = theta1)
arma21 <- arima.sim(model=model,n = 100,innov=rnorm(100,0,sd=sigma))
matplot(arma21,lwd=2,type="l",col="blue",main="Simulated Time Series")

Произведем оценку модели методом максимального правдоподобия

library(TSA)
## Loading required package: leaps
## Loading required package: locfit
## locfit 1.5-9.1    2013-03-22
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-9. For overview type 'help("mgcv-package")'.
## Loading required package: tseries
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
data(arma21)
## Warning in data(arma21): data set 'arma21' not found
mod21 <- arima(arma21,order = c(2,0,1))
mod21$coef
##        ar1        ar2        ma1  intercept 
##  0.4544221 -0.5834244  0.1286499  0.1018173
mod21$sigma2
## [1] 14.78368

Построим прогноз и 95% доверительные интервалы

plot(mod21,n.ahead=12,type='b',xlab='Time',col="blue",lty=3,n1=length(arma21)-20,lwd=1)
abline(h=coef(mod21)[names(coef(mod21))=='intercept'])

Значения прогнозов \(\hat y_t(l)\) и стандартное отклонение прогнозов \(\sqrt{D[e_{t+l}]}\) можно вычислить для все \(l\) c помошью метода predict(). Доверительные интервала любого заданного уровня \((1-\alpha)100\)% можно вычислить , используя квантиль нормального распределения и стандартное отклонение прогнозов \(\sqrt{D[e_{t+l}]}\)

forec<-predict(mod21,n.ahead = 10)
forecpl <- forec$pred+forec$se
forecmn <- forec$pred-forec$se
ff <-cbind(forec$pred,forecpl,forecmn)
matplot(ff,lty= 1,type = "l",lwd = 2,main = "Forecasts +- std value" )

При прогнозировании нестационарного ряда, который сводится к стационарному операцией дифференцирования (метод diff(.)), прогнозируйте сначала стационарный ряд, затем применяйте к прогнозам обратный к дифференцированию оператор суммирования (метод diffinv(.))

Пример прогнозирования смоделированного ряда

Следующий фрейм есть пример прогнозирования смоделированного ряда. Попробуйте разную длину прогноза по истории ряда.

Пример прогнозирования импорта Росии

В прошлой лекции мы построили сезонную ARMA модель для импорта России. Здесь мы посмотрим на его прогнозирование. Сначала повторим построение лучшей модели из прошлой лекции без лишних комментариев. И к ряду из разностей и сезонных разностей, взятых с мая 2006 года по настоящее время подгоним сезонную SARMA(0,0,3)(0,0,1) модель.

y <- diff(importFact)
z <- diff(y,lag = 12)
zz <- z[100:216]
zz <- ts(zz,start= c(2006,5),frequency = 12 )
plot(zz ,type = "l", col = "blue",lwd = 2,main = "Import (mlrd dollar).Russia.Seasonal Differences")

zz <- zz[1:116]
modelimport2 <- arima(zz,order = c(0,0,3),seasonal = list(order= c(0,0,1),period = 12 ),method = "ML")

Как и в прошлом примере построим прогноз на 12 месяцев вперед. Добавим к нему прошлые значения и доверительные интервалы уровня 0.8

getVector <- function(vv,bInd,l)
{
  c <- vector()
  if (bInd > 1)
  {
    c <- rep(NA,bInd-1)
  }
  res <- c(c,vv)
  cc <- vector()
  cc <- rep(NA,l-length(res))
  res <- c(res,cc)
  return(res)
  
}  
lzz <- length(zz)
nforec <- 12
ll <- lzz+nforec
forecImport<-predict(modelimport2,n.ahead = nforec)

vec1 <- getVector(zz,1,ll)
vec2 <- getVector(forecImport$pred,lzz+1,ll)
vec3 <- getVector(zz,1,ll)
q<- qnorm(0.75,0,1)
dPl <- forecImport$pred + forecImport$se*q
dMn <- forecImport$pred - forecImport$se*q
vec4 <- getVector(dPl,lzz+1,ll)
vec5 <- getVector(dMn,lzz+1,ll)
res <- cbind(vec1,vec2,vec3,vec4,vec5)
matplot(res ,  type = c("b","b","b","b","b"),pch = 16,lty=1, lwd= c(2,3,2,1,1),ylab = "Values",xlab="Time",col = c("green","blue","magenta","black","black"))
abline(h=coef(modelimport2)[names(coef(modelimport2))=='intercept'])
      legend("topright",c("Future Values","Forecasts","Time series","Conf.Level"),bty="n",lwd = 2,col = c("green","blue","magenta","black"))

Мы получили прогноз на два года ряда из разностей. Вернемся к исходным значениям взяв обратный к сезонным разностям оператор diffinv. В качестве начальных условий возьмем первые 12 значений из ряда \[y_t=import_t-import_{t-1}\]. Проинтегрируем также доверительные интервалы и такими же начальными условиями.

firsty <- y[1:12]
iiy <- c(z,forecImport$pred)
iPl <- c(z,dPl)
iMn <- c(z,dMn)
iforec <-  diffinv(iiy,lag= 12,differences = 1,firsty)      
iforecPl <-  diffinv(iPl,lag= 12,differences = 1,firsty)      
iforecMn <-  diffinv(iMn,lag= 12,differences = 1,firsty)

Теперь проинтегрируем уже несезонные разности с начальным значением равным первому значению ряда импорта.

  firstImport <-importFact[1] 
  imforec <- diffinv(iforec,lag= 1,differences = 1,firstImport) 
  imforecPl <- diffinv(iforecPl,lag= 1,differences = 1,firstImport) 
  imforecMn <- diffinv(iforecMn,lag= 1,differences = 1,firstImport)
  imforec <- ts(imforec,start = c(1997,1),frequency = 12)
  imforecPl <- ts(imforecPl,start = c(1997,1),frequency = 12)
  imforecMn <- ts(imforecMn,start = c(1997,1),frequency = 12)
  imforecPl[1:length(importFact)]= NA
  imforecMn[1:length(importFact)]= NA
  
  
  plot(imforec,col = "blue",lwd = 2,type = "l", main = "Import Russia Forecasts and Conf.Interval" )
  lines(importFact,lwd = 2,col = "red")
  lines(imforecPl,lwd = 2,lty = 4,col = "black")
  lines(imforecMn,lwd = 2,lty = 4,col = "black")

Увидеть результаты прогнозирования для произвольного момента времени можно в следующем фрейму